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Abstract —We develop a new compressive sensing (CS) inver¬ 
sion algorithm by utilizing the Gaussian mixture model (GMM). 
While the compressive sensing is performed globally on the 
entire image as implemented in our lensless camera, a low- 
rank GMM is imposed on the local image patches. This low- 
rank GMM is derived via eigenvalue thresholding of the GMM 
trained on the projection of the measurement data, thus learned 
in situ. The GMM and the projection of the measurement data 
are updated iteratively during the reconstruction. Our GMM 
algorithm degrades to the piecewise linear estimator (PLE) if 
each patch is represented by a single Gaussian model. Inspired 
by this, a low-rank PLE algorithm is also developed for CS 
inversion, constituting an additional contribution of this paper. 
Extensive results on both simulation data and real data captured 
by the lensless camera demonstrate the efficacy of the proposed 
algorithm. Furthermore, we compare the CS reconstruction re¬ 
sults using our algorithm with the JPEG compression. Simulation 
results demonstrate that when limited bandwidth is available 
(a small number of measurements), our algorithm can achieve 
comparable results as JPEG. 

Index Terms —Compressive sensing, Gaussian mixture models, 
dictionary learning, sparse representation, lensless camera, low- 
rank. 

1. Introduction 

Compressive sensing 0-0 (CS) has led to real applica¬ 
tions, including the single-pixel camera Q, the lensless cam¬ 
era 0 0 video compressive sensing 0 -m depth com¬ 
pressive sensing 0 m hyperspectral compressive imag¬ 
ing GDI T5| , polarization compressive sensing terahertz 
imaging |17| , and millimeter wave imaging |T^ . In this paper, 
we focus on the two-dimensional (2D) image case, though 
similar technique can be used for videos |T9] and other 
band widths. Specifically, we develop our algorithm under the 
lensless compressive imaging architecture Q, p0| , which has 
provided excellent reconstruction images from the compressive 
measurements using simple and off-the-shelf hardware pQ| . 

Diverse algorithms ED -|[26| have been developed for com¬ 
pressive sensing recovery, which plays a pivot role in CS, to 
reconstruct the desired signal from compressive measurements. 
Sparsity, the key ingredient of CS, has been investigated ex¬ 
tensively in these algorithms. The wavelet transformation ED 
p7|-p9| is generally used since it provides the sparse rep¬ 
resentation of an image with fast transformations (thus very 
efficient). A parallel research is using the total variation (TV) 
for CS recovery | [23| , p0| , which provides good results for 
piecewise smooth signals. The aforementioned algorithms do 
not increase the unknown parameters significantly during the 
reconstruction, as usually the wavelet coefficients will have 


similar (or the same) number of the image pixels. Recently, 
researchers have found that by exploiting the sparsity of local 
patches ED’E3’ better results can be achieved. In summary, 
CS recovery algorithms fall into the following three categories: 
1) global basis based algorithm, e.g., using the wavelet trans¬ 
formation, 2) TV based algorithm, and 3) local basis based 
algorithm, i.e., DCT or dictionary learning or denoising based 
algorithms. State-of-the-art CS inversion results have been 
obtained in p4| , which generally lie within the third 
category. In this paper, we propose an alternative inversion 
algorithm that exploits the low-rank property of image patches. 

Generally, the CS recovery is an iterative process in which 
two steps are performed at each iteration i) projecting the 
measurements to the image level, which can be done by the 
majorization-minimization (MM) approach p2| , pS] , or the 
Euclidean projection p6| , or the alternating direction method 
of multipliers (ADMM) p7| , ii) denoising this projected 
image and updating the recovery of the desired image. These 
two steps are performed iteratively until some criterion is 
satisfied. A general framework is developed in ED under 
the approximate message passing (AMP) framework, in which 
diverse denoising algorithms p8| can be plugged in. One 
key difference of the denoising based CS inversion algorithms 
compared with the wavelet based CS inversion algorithms is 
that the former exploits the local sparsity based on (overlap¬ 
ping) patches, as state-of-the-art denoising algorithm is using 
sparse representation of local patches, e.g., p9| . The number 
of coefficients for the patches (under some basis or dictionary) 
is usually larger than the image pixel number (because the 
overlapping patches are used). 

A. Contributions 

The proposed algorithm in this paper also lies in the third 
category mentioned earlier, which is based on the local basis 
(for image patches). Specific contributions of this work can be 
summarized as below: 

• We investigate the low-rank property of image patches 
under the Gaussian mixture model (GMM) framework. 
Different from the low-rank model investigated in 
which requires patch clustering (block matching) as an 
additional step, in our algorithm, each patch is modeled 
by a GMM with different weights corresponding to dif¬ 
ferent Gaussian components, which can be seen as a soft 
clustering approach based on these weights. Furthermore, 
these weights are updated in each iteration. 
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• We develop an general framework using ADMM to 
explore the sparsity (and low-rank property) of patches 
in order to recover the desired signal. 

• If each patch is modeled by a single Gaussian component, 
our GMM degrades to the piecewise linear estimator 
(PLE) | [40| , | [4T| . Therefore, a low-rank PLE algorithm 
is also proposed for CS recovery. 

• We conduct experiments using our proposed algorithm 
and other leading algorithms on the real data captured 
by our lensless camera. This verifies real applications of 
each algorithm. 


B. Organization of This Paper 


We start with the derivation of an ADMM formulation for 
CS recovery to investigate the sparsity of local overlapping 
patches in Section |I^ The proposed low-rank GMM algo¬ 
rithm is developed in Section O and the joint reconstruction 
algorithm is summarized in Section |IV| Extensive results on 
both simulation and real data are reported in Sections [V|[Vl 
Section |VII| concludes the paper. 


II. CS Inversion via Exploiting Sparsity of Patches 


formulation in Section llv] Introducing Lagrange multipliers 
p and parameter p in (jSjTresults in the objective function 

L{x,z,X,p) = ^\\y - Ax\\l + \\\z\\i+ 

i 

+ |^||RiX-B^i||2 (7) 

i 

Define u = {1 /p) p, 

L{x,z,u,X) = i||t/- Aa;||| + All^lli 

+ ^ llR-ja: — + u \\2 + const(8) 

i 

The ADMM cyclically solves the following 3 sub-problems: 

a;*+^ ;= arginin(i||y - Ax\\l + ^ ||Ri® - 'Bz\ + u^\\l) 

X z z 

(9) 

z*+^ := argniin(A||2||i + ^ V ||R.®*+i - Bz, + u%) 

( 10 ) 

;= M* +y^(Ri®*+^ - B 2 ;L^) ( 11 ) 


Under the CS framework, the problem we are solving can 
be formulated as: 

min i||y-Aa;||^+A|| 2 ||i, (1) 

s.t. X = B 2 :, (2) 

where A G is the sensing matrix, x G is the 

desired signal, z is the coefficients which are sparse under the 
basis B. Erom 2 ;, we can recover x easily via B, which can 
be known a priori, (e.g., a wavelet or DCT basis) or learned 
based on x during the reconstruction. 

Considering the image case investigated here, let x denote 
the vectorized image and the sparse representation, 2 :, is now 
modeled on image patches. Therefore, © can be reformulated 
as: 

Rcc = B 2 :, (3) 


where t denotes the iteration index. 

Equation is a quadratic optimization problem and can 
be simplified to 

(A^A + r? y] R7 Ri)x = A^y + p «*)’ 

i i 

( 12 ) 

which admits the closed-form solution, 

a; = (ATA+7?y]R7R0“MA^y+7?y](R7B2i-R7M‘)]. 

(13) 

However, the dimension of {A^ A-\-p large (the 

pixels of the desired image), requiring a high computational 
workload. Alternatively, since (A^ A + 77 Ri) is invert¬ 

ible, the matrix inversion formula can be used to reduce the 
computational workload. As R^ is used to extract i-th patch 
from an image, ^ diagonal matrix 


where R denotes the patch extraction and vectorization oper¬ 
ation and considering each patch, we have 

KiX = Bzi (4) 

with i indexes the patches. 

The problem can be reformulated as: 

min i||y-Aa;||^+A|| 2 :||i, (5) 

S.t. BiX = Bzi. ( 6 ) 

Note this B is shared for all patches for current discussion, 

and in the following analysis, similar patches can be grouped 
together and each group can have its own B. 

Next, we develop an ADMM p7| formulation of 0 to 
solve the problem, which will also be used in our GMM 


R y^R7Ri = diag(ri,..., tn). (14) 

i 

Each of the diagonal entries corresponds to an image pixel 
location and its value is the number of overlapping patches 
that cover that pixel. Therefore, R“^ = diag(rf^,... ,r^^) 
and 

(A’^A + r/R)-i =r/-iR-i 

- 7?-1R-1a’^(I + Ap-^R-^A'^y^Ari-^R-^. (15) 

However, this is not necessarily easy to calculate though R 
can be pre-computed. AR“^A^ needs to be saved for com¬ 
putation. Alternatively, (HI* can be solved by the conjugate 
gradient algorithm | [42| . 

To mitigate this problem, we apply the ADMM again on 
0 and introduce another auxiliary variable w, leading to the 
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following optimization problem: 

min \\\y - +\\\z\\i + '^'^\\RiW ( 16 ) 

s.t. X = w. (17) 

Following this, 

1 71 ^ 

{x,w,z,a) = argmin-||t/ - Aa ;||2 + 7^ X] 




x,w,z,cx. ^ 


/?, 


+ A||^||i+a {x - w) +-\\x - w \\2 


(18) 


which can be simplified to (by setting v = cx//3y. 

7 
2 


{x,w,z,v) = argmini||j/ - Ax\\l + ^ HR^m; - Bz, 

x,w,z,v ^ ' 


.111 


+ ^INIIi + ^11® “ m + v\\\ + const 


(19) 


The optimization of ( p^ consists of the following iterations: 
1 

X 2 


:=wgm.m.^\\y-Ax\\l + ^\\x-w'^ + V^\\l, (20) 


,*+l 

:= argmin ^ V ||RiW; - BzjWl 

i 




(21) 

.t+1 

:= argininA|| 2 :||i -f y Y] ||Ri'u;*+^ - Bzi\\l, 

(22) 

,*+l 


(23) 


For fixed x*'^^ admits the following closed-form 

solution: 

a;‘+i = (A"^A -f /3I)“^[A"^y -h /3(w* - v*)], (24) 

which can be simplified to 

a:*+i = - /?-iA'^(I -f A/3-^A'^)-^A/3-^) 


X [A'^y + l3{w* -■n*)], 


(25) 


For the case considered in our work (as implemented in the 
lensless camera j^), A is the permuted Hadamard matrix and 
thus AA^ is an identity matrix: 


X 


^+1 




A^A 


/3 -f 1 
= (m* - -w* 


{P + 1)^ 

+ {w^ - v^) - 


[A^y + 


A^y ^ ^ A^ A{w^ — v^) 


/3 + 1 


A~^ [y — A(w^ — v^)) 

^+1 ■ 


(26) 


Similarly, for fixed admits the following 

closed-form solution: 


w 


t+i 


= U^R7Ri+/3I 


/3(x*+^+v*)+7)J2R7Bzj 


(27) 


Recall that ^ Ri is a diagonal matrix R = 
diag(ri,..., tat), thus w can be computed element wise via 


(+1 _ 

“ vr„ + /3 


(28) 


where [-J^ denotes the n-th entry of the vector inside [ ]. 

Similar to (10), (22) can be considered as a dictionary 
learning model, where B is the dictionary. If the orthonormal 
transformation is used (e.g., the DCT), z can be solved by the 
shrinkage thresholding operation p^, p7|. 


Since the key of this algorithm is to investigate the sparsity 
of the local overlapping patches, we term this framework as 
SLOPE (Shrinkage of Local Overlapping Patches Estimator), 
where the ‘local’ stands for the local basis rather than the 
global basis such as wavelet. The ADMM-SLOPE is summa¬ 
rized in Algorithm 


While good results have been obtained using similar ap¬ 
proaches ED, as in Algorithm in the next section, we 
develop a low-rank GMM framework imposed on the patches 
and a full formulation of the proposed algorithm is presented 
in Section lEl 


Algorithm 1 ADMM-SLOPE 

Require: Measurements y, sensing matrix A, r], A}. 
1: Initial x^w^v to all 0. 

2: for t = 1 to Max-Iter do 
3: Update X by Eq. ( [26|) . 

4: Update w by Eq. \21) . 

5: Update 2 : by shrinkage operator. 

6: Update V by Eq. ( [2^ . 

7: end for 


III. The Gaussian Mixture Model 

The Gaussian mixture model (GMM) has been re¬ 
recognized as an advanced dictionary learning approach and 
has achieved excellent results in image processing ED 
and video compressive sensing p^ , Recall the image 
patches X G extracted from the 2D image, where the 

patch size is \/P x \/P and there are in total Np patches. For 
i-th patch Xi, it is modeled by a GMM with K Gaussians [ [44| : 

K 

ajj ~ ^ TTkM{p.k, Sfe) (29) 

k=l 

where represent the mean and covariance matrix 

of k-thc Gaussian, and denotes the weights of these 

Gaussian component. 

In this paper, we further impose the GMM is low-rank and 
now the model in becomes ( [^ and the problem to be 
solved becomes 


min -||y-Aa;||^ 

(30) 

K 


s.t. ajj ~ ^ TTkAfijik, Sfe) 

(31) 


k=l 
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where Xi denotes i-th patch from x, which is an vectorized 
image. symbolize the low-rank GMM. 

The following problem is to estimate this low-rank GMM. 
Recalling Section [T| we review that the CS recovery is an 
iterative two-step procedure. In each iteration, one can get 
an estimate from the projection of the measurements (details 


We further define 

= Vfc = l,-' 
Next, we define a new random vector 

K 


,K. 


discussed in Section [IV| ). We hereby learn a (full rank) 
GMM from this estimate and then proposing the eigenvalue 
thresholding approach to derive the low-rank GMM based on 
this full rank GMM. 




(41) 


(42) 


k=l 


which is modeled by a low-rank GMM, parameterized by 


A. Low-Rank GMM 

In the following, we provide a motivation for the next step 
in the algorithm. The random vector Xi in ( [^ (modeled 
as a GMM) can be written as, dropping the subscript i for 
simplicity, 

K 

= '^T^kgk, (32) 

k=l 


B. Update Estimate via the Low-Rank GMM 


X 


9k 


9k — ^kQk + Mfc) 


where G 




= argrnin{-||r - ^ rk\\T\\^}, VA: = 1, • 


^k = U/jA/cU^, 

A/c = [Ai,...,Ap]. 

We impose that 7 /. < P via 

A/c = [Ai,..., At^^ , 0], 

Xi = max(Ai - A^^+i, 0), Vi ^ 


^/c = U/cA/cU;. 


where each is a random vector of multivariate normal 
distribution, given by 


Given an estimated image x, the GMM in ( [2^ can 
be learned via the Expectation-Maximization (EM) algo¬ 
rithm | [43l , | [44| based on overlapping patches. Then for each 
Gaussian component, we adopt the EVT to the covariance 
matrix to obtain the low-rank GMM {/i;., Eollowing 

this, the estimated image x or patches Xi can be updated via 
this low-rank GMM, to x or xp Dropping the subscript i, 
given X, the conditional distribution for x maybe evaluated as 


p{x\x) = 


p{x)p{x\x) 


(33) 


The random vector can be decomposed into independent 
random variables of normal distribution p3l, p31 as follows 


J p{x)p{x\x)dx 
Since x is a low-rank version of x, we assume: 


X = X 


n, 


(43) 


(44) 


(34) 


G 7 /c is the rank of and 


where n ~ A/'(0, E) is modeled as an additive Gaussian noise, 
thus 


is a vector whose jk components are independent random 
variables. 

In order to reduce noise, we use a model in which has 
a small number of independent random components, i.e., we 
require 7 /c to be small. This is equivalent to requiring have 
a reduced rank. More specifically, introducing the parameter 
T/c, we solve the following minimization problem ||46l: 


,K. 

(35) 


Plugging 

pi^x) 


p{x\x) ~ M{x,'Ei) 
into l |43| i, we have 

p{x)p{x\x) 
f p{x)p{x\x)dx 

^ J2k=l T^kJ^iilk^ Sfc) X J\f{x, E) 

/ S;) X V(i, E)d* 

K 

fe=i 


where || • ||p is the Erobenious norm, and || • ||* is the nuclear 
norm (sum of the singular values). It is shown in |[46| that the 


solution to (35) can be readily obtained by a shrinkage on the 
singular values (which are similar to the eigenvalues) of 
Specifically, consider 


which is an analytical solution with | [44| 


(36) 

(37) 


(38) 

1,...,P. (39) 


<t>k 

flk 

^k 


'7^/-A/'(x; E + ^k) 

(E 

^k 


-1 


^/c(E + E/c) 


= ri/e(E fl^) 

— ^/c(E + '^k) ~ P'k) “b A/c- 


(45) 

(46) 

(47) 

(48) 

(49) 

(50) 

(51) 


And we term this as the eigenvalue thresholding (EVT). 
Eollowing this, E/c is obtained by 


Note that is low-rank obtained via EVT from but 
by adding E (= cr^Ip), (E + E/^) is invertible. While (48) 
provides a posterior distribution for x, we obtain the point 
estimate of x via the posterior mean: 


K 


(40) 


E[a;] = (pkt^k 


(52) 


k=l 
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which is a closed-form solution. 

The procedure of learning and updating the GMM can be 
summarized as below: 

• Step 1: Lean a GMM (not low-rank) {tt/c, '^k}k=i 
EM from an estimate of x, which can be obtained from 
1ST, GAP or ADMM described below in Section |IV| 

• Step 2: For each Gaussian component, derive the low- 

rank version eigenvalue value thresh¬ 

olding via 

• Step 3: Update the estimate of the image by x using 

(148](-(152](. 


C Degrade to the Piecewise Linear Estimator 

The piecewise linear estimator (PLE) proposed in ED 
has demonstrated excellent performance on diverse image 
processing tasks. If each patch is considered drawn from a 
single Gaussian distribution, our GMM degrades to the PLE 
and the weights tt/c (or (j)k) are not required. Furthermore, 
the update equation of x will become the Winer filter. The 
MAP-EM procedure proposed in ED can still be used to 
determine which Gaussian each patch lies in and to estimate 
the denoising version of the each patch. However, this method 
has been shown that it is very sensitive to the initialization and 
selecting K is critical to the performance of the method. We 
compare our proposed algorithm with PLE by experiments in 
Section IV-DI 

It is worthing nothing that, even using PLE, the eigenvalue 
thresholding method used to obtain the low-rank Gaussian 
model is first proposed in this paper. In this case, the PLE 
model is very similar to the NLR-CS p^ , where the low- 
rank is imposed on each cluster of patches, while in the PLE, 
the low-rank is imposed on the patches belonging to the same 
Gaussian component; this can also be seen as a cluster. 

Next, we review the MAP-EM algorithm proposed in ED 
and adopt it to the current context. In the E-step, assum¬ 
ing that the estimates of the low-rank Gaussian parameters 
^k}k=i known, (following the previous M-step), for 
each patch, one calculates the MAP estimates 0^ of all the 
Gaussian models and selects the best Gaussian model ki to 
obtain the estimate of the patch Xi = 6\\ In the M-step, 
assuming that the Gaussian models selection ki and the signal 
estimate are known (following the previous E-step), one 

updates the Gaussian models then impose 

them to be low-rank. 

• E-step: Signal Estimation and Model Selection: 

For each image patch i, the signal estimation and the 
model selection are calculated to maximize the log a 
posteriori probability \ogp{xi\xi)\ 


{xi,ki) = argmaxlogp(0|®i,/ife,i:fc) (53) 

d,k 

Recall that we consider x is a low-rank version of Xi and 
Xi = Xi-fn, n ^ Af{0,a‘^lp) (54) 


Therefore: 


(xi.ki) = argmax (logp(cci|0,cr^lp) 

+ logp(0|/ifc,Sfc)) (55) 

= argmin (a~‘^\\xi - 6\\‘^ + 0.5log |Sfe| 

G,k \ 

+i^ - - P-k)) ( 56 ) 

This maximization is first calculated over 6 and then over 
k. Given a prior Gaussian signal model 6 r\j 
9 can be estimated by the posteriori mean 

0*= = tk{tk+aHp)-^Xi (57) 

The best Gaussian model ki that generates the maximum 
MAP probability among all the models is then selected 
with the estimated x^ 

ki = argmin (cr'^llccj - + 0.51og |Sfe| 

+{e>i - - ii,)'y5S) 

The signal estimate is obtained by plugging in the best 
model ki in the MAP estimate 

Xi = df- (59) 


M-step: Model Estimation: 

In the M-step, the Gaussian model selection ki and the 
signal estimate Xi of all the patches are assumed to be 
know (derived from the 1ST, GAP or ADMM as shown in 
Section W). The parameters of each Gaussian model are 
estimated with the maximum-likelihood (ML) estimate 
using all the patches in the same Gaussian model: 


= arg max logp 1^*., f2fc) (60) 


where Ck denotes the ensemble of the patch indices i that 
are assigned to the k-\h Gaussian model and 




(61) 




^k — P'k) • 




Return to the low-rank model proposed in this paper. The E- 
step is same as above and one more step is added in the M- 
step. The full rank (not low-rank) Gaussian models are first 
estimated via ([^-([^, and then each Gaussian model is im¬ 
posed to be low-rank by thresholding the eigenvalues via ( [3^ - 
( [40| ). The new low-rank PLE algorithm can be summarized 
into the following 3 steps: 


Step 1: Signal estimation and model selection by ED- 


Step 2: Model update for the (non low-rank) Gaussian 
models by 

Step 3: Estimate the low-rank Gaussian models 


{P'k^'^k}k=i from {tik,'^k}k=i via M-(41). 
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IV. The Joint Reconstruction Algorithm 

Sectioning presents an algorithm to obtain a better estimate 
(or a denoised version) of the signal x given an initial estimate 
utilizing the low-rank GMM. In this section, the GMM will be 
wrapped into our joint reconstruction algorithm by different 
update methods to get the initial estimate, which can be 
considered as projecting the measurement y to the image plane 
X. This is obtained by minimizing the following objection 
function 

= ||y - Ax||^. (63) 


Diverse algorithms have been proposed and we review two 
of them below and develop an ADMM formulation in Sec¬ 


tion Other approaches, for example, the TwIST p7| can 
also be used. 


A. Iterative Shrinkage Thresholding 

By using the majorization-minimization approach to 
minimize J{x), we can avoid solving a system of linear 
equations. At each iteration t of the MM approach, we should 
find a function Gt{x) that coincides with J{x) at x^ but 
otherwise upper-bounds J{x). We should choose a majorizer 
Gt{x) which can be minimized more easily (without having 
to solve a system of equations). The Gt{x) is defined as 

Gt{x) = \\x — Ax\\l + {x — x^)^ {(1 — A){x — x^), (64) 

where I denotes the identity matrix and ( must be chosen to 
be equal to or greater than the maximum eigenvalue of A. 
For the Hadamard sensing matrix used in our camera, the 
maximum eigenvalue of A^A is easily obtained. The update 
equation of x^ in this Iterative Shrinkage Thresholding (1ST) 
algorithm p5| is given by: 

^t+i = a;* + 1(y - Aa;*). (65) 

B. Generalized Alternating Projection 

The Generalized Alternating Projection (GAP) algorithm 
proposed in pb} , which enjoys the anytime property and has 
been demonstrated high performance in video compressive 
sensing has the following update equation by using the 
Euclidean projection: 

=a:* + A'^(AA"^)-i(y-Aa;*). (66) 

Under some condition of the sensing matrix A, as the 
Hadamard matrix used in our system, AA^ is the identity 
matrix and thus ( [66| ) is same as ( [65] ) with C = 1- 

In addition to ( [6^ , aiming to speed-up the convergence, the 
authors in | [36| have proposed the accelerated update equations 

= cc^-hA“^(AA'^)-^(y^-Ax^), (67) 

y* = y*-i + (y - Aa;*"^). (68) 

Better results have been achieved in our experiments using 
this accelerated GAP. 


Algorithm 2 LR-GMM-SLOPE 

Require: Measurements y, sensing matrix A. 

1: Initial x. 

2: for f = 1 to Max-Iter do 

3: Update a; by 1ST or GAP ^ or ADMM (|7T]i. 

4: Update related parameters in 1ST, GAP or ADMM. 

5: Learn a GMM (not low-rank) from x. 

6: Obtain the low-rank GMM via eigenvalue shrinkage 

thresholding ( [^ . 

7: Update X by the low-rank GMM using expectation 

in ( [5^ . 

8 : end for 


Algorithm 3 LR-PLE-SLOPE 


Require: Measurements y, sensing matrix A. 

1: Initial x. 

2: for f = 1 to Max-Iter do 

3: Update X by 1ST or GAP ^ or ADMM 

4: Update related parameters in 1ST, GAP or ADMM. 

5: Update the Gaussian models (not low-rank) from x via 


6: Obtain the low-rank Gaussian models via 

7: Update X by the low-rank Gaussian models using ([57 )- 

•HU. 

8 : end for 


C. An ADMM Formulation 


Under the GMM framework, we don’t have the sparse 
variable 2 ; as in 0, the objective function can be formulated 
as: 

X = argmin -||y - Aa :||2 + iY] (^9) 

i 

where x is obtained by the low-rank GMM model. Following 
the procedure in ( p^ by introducing the auxiliary variable 
{tu, v}, we have 

{x,w,v) = argmini||y - Ax\\l + | ||RjM; - Xi\\l 

x,w,v ^ ^ 

I 

+ - tu + 1;||2 + const (70) 

The optimization of ( |70| ) consists of the following iterations: 
:= argmin Ly - Aa ;||2 + ^||a; - + r;*|| 2 , (71) 

X 2 2 

f + (72) 

^t+1 _ ^^t+1 _ ^t+ly (-73^ 


where the update of x is given by the low-rank GMM in ([^ 


(52). Under the sensing matrix considered here in our work, 
AA^ = I, the solution of (71_) is given by (^). Eq. (72) can 
be solved by 




( 74 ) 
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Similar to (28), can be solved element-wise but in one 

shot. 

The proposed low-rank GMM algorithm, integrated with the 
three approaches to update x (projecting y to x), constitutes 
the LR-GMM-SLOPE algorithm summarized in Algorithm 
Similarly, when the low-rank constraint is imposed on the 
PLE, we obtain the LR-PLE-SLOPE in Algorithmic 


Eollowing the formulation of the lensless camera 0 , the 
permuted Hadamard matrix is used as the sensing matrix. Each 
image is resized to 256 x 256 and these images are shown in 
Eigure[C The CSr is defined as: 


CSr = 


number of rows in A 
number of columns in A 


(75) 


V. Simulation Results 


We test the proposed algorithm on simulation datasets 
with 2D images. The proposed algorithm is compared with 
other leading algorithms 1) TVAL3 2) GAP based on 
wavelet p 6 |, 3) DAMP with BM3D denoising, and 4) 
NLR-CS ||33|, which explores the low rank of similar patches. 
State-of-the-art results have been obtained by j^, p4| . The 
Gaussian components in our mixture model is set to AT = 6 for 
all the experiments and the analysis of this number is provided 


in Section |V-C| When updating x, accelerated GAP in ( [67] ) 
is used and the comparison of different approaches is shown 


in Section V-B We obtained the low-rank GMM by setting 
the rank of each Gaussian component learned via the EM 
algorithm to the half of the full rank 7 /. = 0.5P, where P 
= 64 for the patch size 8 x 8 used in this paper. The proposed 
algorithm is further compared with JPEG compression in 
Section IV-AI 



Fig. 1. Images used for CS experiment, {barbara, boat, cameraman, foreman, 
house, lena, monarch, parrot}. 


Proposed, psnr: 26.039 



TVAL3, psnr: 23.2817 



DAMP, psnr: 24.8019 




GAP, psnr: 23.3256 



NLR-CS, psnr: 25.7166 



ground truth 


Fig. 2. Reconstruction results of different algorithms at CSr= 0.1, image 
size 256 x 256. 


where the number of columns in A is equivalent to the total 
pixel number of the image. The sensing matrix is constructed 
from rows of a Hadamard matrix of order N = 2^^. The 
columns of the Hadamard matrix is permuted according to a 
predetermined random permutation (the same permutation is 
used in the real data captured by our lensless camera). Eor each 
CSr, we use the first CSrxA rows of the column-permuted 
Hadamard matrix as the sensing matrix. By selecting some 
other rows of the permutated Hadamard matrix, we can get 
better results | [47| . However, here we just select the top rows 
to be consistent to the implementation of our lensless camera. 
Note that in this case, AA^ is an identity matrix and it is very 
fast to use accelerated GAP update for x in •Hz). We observed 
that best results are obtained by the LR-GMM-SLOPE with 
GAP updates of x and these results are reported in this section. 
Eor the comparison of 1ST, GAP and ADMM, please refer to 
Section I V-B I 

Since when CSr=0.1, good results have been achieved for 
most images (see Eigure for one example), we here spend 
more efforts on the extremely low CSr, in particular CSr< 0.1, 
which may be of interest to the real applications that are used 
to detect anomalous events | [48| without caring too much about 
the image quality. The results are summarized in Table We 
can observe that best results are obtained by the proposed 
aglorithm, NLR-CS or D-AMP. When CSr is less than 0.1, 
the proposed algorithm usually provides best results except 
“foreman”, where NLR-CS is the best. We also observed that 
NLR-CS is very sensitive to the parameters and sometimes the 
PSNRs are not linearly increasing as the CSr increases, while 
the other algorithms including the proposed do not have this 
problem. On average, our proposed algorithm works best when 
CSr< 0.1. Though did not reported here, when CSr> 0.1, the 
proposed algorithm also provides comparable or better results 
than D-AMP and NLR-CS. But we need to tune the noise 
parameter used in n, while for all the results presented here, 
we set it to the same value E = 10“^I. In addition, we need 
to tune the rank thresholds 7 /^, for which we have observed 
that a higher CSr requires a larger 7 /.. 

In addition to the grayscale images tested above as in other 
papers, we also conduct our proposed algorithm on the RGB 
images with results shown in Table |I^ Three sensors are 
simulated to capture the R, G and B components of the image. 
The “book” scene corresponds to the real data captured by 
our lensless camera. Again, our proposed algorithm provides 
the best results. One example with CSr= 0.03 is shown in 
Eigure It can be seen that our algorithms provide more 
details than NLR-CS and D-AMP; both of them presents 
“blob” artifacts. 
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TABLE I 

Reconstruction PSNR (dB) of different images with diverse algorithms at various CSr. 


Image 

Method 

CSr= 0.03 

CSr= 0.04 

CSr= 0.05 

CSr= 0.06 

CSr= 0.07 

CSr= 0.08 

CSr= 0.09 

CSr= 0.1 


Proposed 

21.7416 

22.7603 

23.5292 

24.1613 

24.6844 

25.2192 

25.6663 

26.0390 


NLR-CS 

20.4887 

21.8187 

20.9706 

21.6355 

24.5237 

25.4407 

25.5486 

25.7166 


d-amp 

19.5424 

20.4412 

21.6364 

22.2731 

23.0846 

23.8539 

24.5367 

25.0486 


GAP-w 

20.0870 

20.9902 

21.6160 

22.2509 

22.6795 

23.0811 

23.4164 

23.7371 

barbara 

TVAL3 

18.1065 

19.3655 

20.6883 

21.4713 

22.0738 

22.8414 

23.0771 

23.4208 


Proposed 

22.7729 

22.8600 

24.6319 

25.2106 

25.9662 

26.5287 

27.1382 

27.7893 


NLR-CS 

21.8342 

21.4521 

24.3866 

16.9015 

22.2784 

24.5024 

23.7634 

23.7528 


D-AMP 

19.7960 

20.9017 

21.9684 

22.8501 

23.5793 

24.2540 

24.9395 

25.5490 


GAP-w 

20.6312 

21.2215 

21.8709 

22.3753 

22.8921 

23.3859 

23.7879 

24.1355 

boat 

TVAL3 

18.2669 

19.4750 

20.2451 

21.5041 

21.9831 

22.3542 

22.7917 

23.1918 


Proposed 

21.3706 

22.1542 

22.8186 

23.4861 

24.0118 

24.4425 

24.8875 

25.1473 


NLR-CS 

12.0313 

16.6251 

16.0262 

17.7892 

17.8507 

18.7355 

19.2338 

21.5345 


D-AMP 

18.3514 

19.3660 

20.3618 

21.2030 

22.0050 

22.7786 

23.3635 

24.0456 


GAP-w 

19.1932 

19.7725 

20.2935 

21.1823 

21.5875 

21.8644 

22.1527 

23.8935 

cameraman 

TVAL3 

17.8787 

18.6441 

20.1835 

20.2041 

20.8832 

21.2041 

21.2917 

21.4903 


Proposed 

28.8157 

29.8322 

30.5442 

31.5101 

32.0363 

32.5592 

33.2439 

33.5244 


NLR-CS 

29.7873 

30.4427 

31.9318 

33.4732 

34.0312 

34.8470 

34.9993 

34.4573 


D-AMP 

22.3218 

24.4925 

26.1075 

27.5484 

28.9518 

30.2165 

31.4011 

32.3636 


GAP-w 

24.5972 

25.5398 

26.2436 

26.8443 

27.3006 

27.8175 

28.2215 

26.6864 

foreman 

TVAL3 

19.3428 

20.9050 

23.1337 

24.1152 

24.9492 

25.42780 

25.9380 

26.5072 


Proposed 

26.7712 

28.5080 

29.5543 

30.1185 

31.0571 

31.5114 

32.1238 

32.6598 


NLR-CS 

25.8497 

28.1219 

30.6011 

27.6256 

30.1692 

31.5297 

28.8471 

29.2505 


D-AMP 

21.7965 

23.7841 

25.2089 

26.6130 

27.8643 

28.9586 

30.0785 

31.1040 


GAP-w 

23.0012 

23.8285 

24.5941 

25.3161 

25.8258 

26.3531 

26.8563 

27.2498 

house 

TVAL3 

19.0674 

20.6872 

21.7140 

23.4273 

23.7382 

24.0901 

24.5912 

25.1538 


Proposed 

22.9046 

23.9496 

24.6221 

25.3661 

25.9910 

26.3633 

26.9981 

27.3694 


NLR-CS 

22.6851 

22.6874 

24.8531 

23.2517 

21.2840 

24.2211 

25.7213 

27.2040 


D-AMP 

20.4629 

21.7906 

22.5929 

23.6507 

24.2494 

24.7895 

25.3561 

25.8032 


GAP-w 

20.7381 

21.4840 

22.2093 

22.6754 

23.0767 

23.5415 

23.9218 

24.2559 

lena 

TVAL3 

19.1813 

19.7752 

20.7421 

21.7209 

22.1211 

22.6207 

23.0676 

23.6665 


Proposed 

19.2072 

20.1870 

21.3428 

22.3022 

22.9918 

23.6285 

24.2536 

24.7206 


NLR-CS 

16.1196 

17.5697 

17.1601 

14.6306 

15.5702 

18.8192 

20.8003 

20.8957 


D-AMP 

16.8849 

17.6731 

18.8636 

19.7654 

20.8845 

21.8400 

22.7742 

23.4913 


GAP-w 

17.3904 

18.0572 

18.8062 

19.3374 

19.8313 

20.2766 

20.7840 

21.1890 

monarch 

TVAL3 

16.2031 

17.4510 

18.0521 

18.6358 

19.0864 

19.6194 

19.9829 

20.4556 


Proposed 

23.1402 

24.1143 

24.9261 

25.5033 

26.4840 

26.9776 

27.4935 

28.1257 


NLR-CS 

20.4401 

20.9168 

22.6560 

22.1715 

22.1582 

22.4706 

24.1134 

26.3153 


D-AMP 

19.5304 

20.8630 

21.9746 

22.8995 

23.9191 

24.6991 

25.5513 

26.5082 


GAP-w 

20.9576 

21.9057 

22.5781 

23.1853 

23.7812 

24.2326 

24.7597 

25.1206 

parrot 

TVAL3 

18.2386 

19.7337 

21.4420 

22.1800 

21.8906 

22.6345 

22.9478 

23.5418 


Proposed 

23.3141 

24.4109 

25.2446 

25.9507 

26.6528 

27.1457 

27.6911 

28.1719 


NLR-CS 

21.1545 

22.4543 

23.5732 

22.1848 

23.4832 

25.0708 

25.2534 

26.1408 


D-AMP 

19.8358 

21.1640 

22.3393 

23.3504 

24.3173 

25.1738 

26.0001 

26.7392 


GAP-w 

20.8245 

21.5999 

22.2766 

22.8427 

23.3212 

23.7845 

24.2015 

24.5660 

average 

TVAL3 

18.2857 

19.5046 

20.7751 

21.6573 

22.0907 

22.5990 

22.9610 

23.4285 


TABLE II 

Reconstruction PSNR (dB) of RGB images with diverse algorithms at various CSr. 



Method 

CSr= 0.03 

CSr= 0.04 

CSr= 0.05 

CSr= 0.06 

CSr= 0.07 

CSr= 0.08 

CSr= 0.09 

CSr= 0.1 

Proposed 

22.0258 

23.0449 

23.8714 

24.4931 

25.0737 

25.6280 

26.0476 

26.4061 

NLR-CS 

12.9754 

13.8812 

19.0401 

17.2073 

17.5413 

18.9149 

22.0660 

21.8692 

D-AMP 

17.3405 

18.9956 

20.6016 

21.8303 

22.8833 

23.8037 

24.6248 

25.1828 

GAP-w 

18.9621 

19.8426 

20.5282 

21.0880 

21.5961 

21.9764 

22.2893 

22.5441 

TVAL3 

15.9549 

17.0403 

18.0474 

18.9911 

19.4949 

20.1248 

20.4482 

20.8150 

Proposed 

22.4159 

23.4245 

24.4111 

25.3072 

26.1026 

26.7247 

27.3431 

27.8475 

NLR-CS 

15.5817 

19.0300 

19.9201 

19.4926 

21.1154 

23.7973 

24.8067 

23.6006 

D-AMP 

17.3405 

18.9956 

20.6018 

21.8303 

22.8833 

23.8037 

24.6248 

25.1828 

GAP-w 

18.9621 

19.8426 

20.5282 

21.0880 

21.5961 

21.9764 

22.2893 

22.5441 

TVAL3 

15.9549 

17.0403 

18.0474 

18.9911 

19.4949 

20.1248 

20.4482 

20.8150 


A. Compare to JPEG Compression 

We now compare LR-GMM-SLOPE under the compressive 
sensing framework with the JPEG compression, which is based 
on the sparsity of the DCT coefficients of 8 x 8 blocks (non¬ 
overlapping patches). We first use an PNG file as the ground 
truth and then use the script within MATLAB “imwrite(-)” 
by choosing 8-bits ‘jpeg’ compression with different qualities 
(100 denotes the highest quality). We treat the quality 100 
as the standard full file size. Eor the ‘Barbara’ image we 
used here, PSNR = 58.47dB (w.r.t. the PNG file) and the 
file size is 45.6KB at quality 100. The compressed image is 


obtained by changing the compression quality from 1 to 100 
and we compare the file size with the full size at quality 100, 
computing the CSr used in this paper. 

Table summarizes the results of JPEG compression 
compared with the results obtained by our algorithm. This is a 
rough, high level comparison because JPEG also performs an 
entropy encoding after the DCT transform and quantization, 
while in our method, the number of compressive measure¬ 
ments is compared with the number of total pixels, and we did 
not consider the entropy coding on quantized measurements. 
We intend to take the effect of the entropy encoding in JPEG 
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Proposed, CSr: 0.03 



NLR-CS 



D-AMP 



GAP wavelet 



TVAL3 




ground truth 


Fig. 3. Reconstruction results of different algorithms at CSr= 0.03. 


TABLE III 

JPEG COMPRESSION AT DIFFERENT QUALITIES COMPARED WITH THE 
PROPOSED COMPRESSIVE SENSING RECOVERY 


JPEG Compression 

CS Reconstruction 

Difference 

Quality 

Size 

(bytes) 

PSNR 

(dB) 

CSr 

PSNR 

(dB) 

PSNR 

(dB) 

1 

1,563 

22.0683 

0.0334 

22.1507 

-0.0823 

3 

1,716 

22.7278 

0.0367 

22.4125 

0.3154 

5 

2,153 

24.8416 

0.0460 

23.2433 

1.5983 

7 

2,559 

26.0924 

0.0547 

23.8426 

2.2548 

9 

2,946 

26.9623 

0.0630 

24.3308 

2.6315 

10 

3,141 

27.2853 

0.0672 

24.5945 

2.6908 

12 

3,494 

27.8379 

0.0747 

24.9514 

2.8866 

14 

3,815 

28.2981 

0.0816 

25.3233 

2.9748 

16 

4,160 

28.6912 

0.0890 

25.6880 

3.0032 

18 

4,478 

29.0306 

0.0958 

25.9230 

3.1077 

20 

4,752 

29.3411 

0.1016 

26.1267 

3.2144 



JPEG compression, quality: 1, PSNR: 22.0683 


Proposed, CSr: 0.033427, PSNR: 22.1507 




JPEG compression, quality: 20, PSNR: 29.3411 Proposed, CSr: 0.10163, PSNR: 26.1267 


Fig. 4. Example images: JPEG compared with the proposed algorithm. 




out of the comparison by computing the JPEG compression 
ratio as compared to the quality 100. When the compression is 
high (lower CSr), the gap between our approach and the JPEG 
compression is small. It is worth noting that when CSr=0.0334, 
our algorithm performs better than JPEG. When the compres¬ 
sion gets lower, the gap becomes larger. One possible reason is 
that when JPEG is performed on the image, the ground truth is 
available and it is very easy to capture useful information from 
the truth. However, under the compressive sensing framework 
and using the current algorithm, increasing a few number 
of measurements can help the reconstruction, but not that 
significantly. Example images can be found in Eigure|^ It can 
be seen that JPEG compression has obvious block artifacts 
while the results of the proposed algorithm become better 
progressively with increasing number of measurements. 

Eurthermore, in JPEG compression, if we lose some bits, 
we may not be able to decode entire blocks. By contrast, 
in our compressive sensing framework, if we lose some 
measurements, we can still reconstruct the image, maybe not 
at a high fidelity. 


Parrot 

40 I-^- I - 1-^- 



15'-'-'-'-'-'-'-'-'-'- 

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 


CSr 

Fig. 5. One example (parrot) to compare different updating rules (1ST, Acc- 
GAP, ADMM) of £c. 


B. Comparison of Different Update Rules for x 


We provide three approaches in Section |IV-A to update x 
in order to minimize the objective function in (|63|). Now we 
compare these three updates via experiments. We emphasize 
again that we are using the permuted Hadamard matrix as the 
sensing matrix and thus AA^ is an identity matrix. Therefore, 
updating x via GAP in ( [^ is same as updating x via 
1ST in ( [65] ). However, the accelerated GAP in ( [67] ) provides 
best results in our experiments. Without tunning the ADMM 
parameters carefully, we compare these three update methods 
with different images at various CSr, and one example is 
shown in Eigure It can be observed that the accelerated 
GAP update always provides the best result and when CSr is 
low, 1ST is better than ADMM. When CSr is getting larger, 
ADMM becomes better than 1ST. Because of this, all the 
results reported in this paper is generated by the accelerated 
GAP update. 


C. Different Number of Gaussian Mixture Components 

One problem of using GMM is how to set the component 
number K. As we are using the mixture model, each patch 
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Fig. 6. Reconstruction PSNR with different number of GMM components 
{K). Barbara is used with CSr = 0.1. 


is represented by the posterior distribution, another GMM. 
Therefore, selecting this K is not as critical as in the PLE ED- 
An alternative way to infer this K is utilizing the manifold 
factor analysis model as developed in | [44| . Hereby, we in¬ 
vestigate this point empirically by using the “barbaba” image 
as used before with different number of iT G [2,20]. The 
results at CSr= 0.1 are shown in Figure It can be seen that 
our algorithm is not sensitive to this AT, since the standard 
deviation of the PSNRs with different K is only 0.04265dB 
compared with the mean value 26.0707dB. 


Barbara 



0.05 0.1 0.15 0.2 0.25 


CSr 

Fig. 7. One example (barbara) to compare LR-PLE-SLOPE (20 Gaussian 
models) with LR-GMM-SLOPE {K = 6). 


D. GMM PLE and Computational Time 

As mentioned earlier, when we consider that each patch 
is drawn from a single Gaussian component, the proposed 
approach degrades to the PLE. We verify the performance 
of LR-PLE-SLOPE compared with the LR-GMM-SLOPE in 
Figure [7] It can be seen that the GMM always performs better 
than the PLE at lower CSr. When CSr is getting larger, they 
start to perform similarly. 

Regrading the computational time, our algorithm is similar 
to NLR-CS. If a warm start is used to initialize the x, we 
can obtain a good reconstruction within 20 iterations. One 
256 X 256 grayscale image reconstruction at CSr= 0.1 takes 
around 1 minute on an i7CPU with 24G RAM. Similar time is 
required for the LR-PLE-SLOPE but it needs more memory. 
While the most time consumption of LR-GMM-SLOPE is 
the EM training of GMM, the LR-PLE-SLOPE requires a 
long time for the model selection, and it usually needs more 


Gaussian components (i.e., 20) than the GMM to get good 
results. 

VI. Real Data Results for the Lensless Camera 

We now verify our proposed algorithm on the real data 
captured by our lensless camera Q, which is composed of 
an aperture assembly and a single sensor (a photodiode) to 
capture grayscale images; it can also be a RGB sensor to 
capture color images. The aperture assembly implements the 
sensing matrix and we programmed it to be the permuted 
Hadamard matrix. By capturing the scene with different 
sensing matrices, we obtain the measurement vector y. We 
implemented the aperture assembly with a transparent LCD 
and thus we can control the image resolution by merging the 
neighboring pixels. 

A. Gray-Scale Images 

We first consider the case with gray scale sensor and the 
image resolution of 128 x 128. To capture compressive mea¬ 
surements, we use a sensing matrix which is constructed from 
rows of a Hadamard matrix of order = 2^^. Each row of the 
Hadamard matrix is permuted according to a predetermined 
random permutation. The scene is composed of a photo printed 
on a paper and we capture the measurements of this photo. 
Example results using different numbers of measurement are 
shown in Fig. We also compare the five algorithms used in 
the simulation. It can be seen that, similar to the simulation, 
our proposed algorithm provides best result when CSr is 
small. Especially, at CSr = 0.05 and 0.1, our algorithm can 
present many details of the face, for example, the left eye 
of “Lena”. D-AMP introduces some “blob” noise because 
the BM3D denoising approach is used. Though NLR-CS can 
provide good results at CSr = 0.05 and 0.1, it introduces some 
unpleasant artifacts when CSr is from 0.15 to 0.25. We also 
tried the algorithm (sHM) proposed in [ [4^ , where a Bayesian 
model is developed to investigate the tree-structure in wavelet. 
Surprisingly, sHM now works better than TVAL3 and GAP, 
mainly due to the following two reasons. Firstly, the tree 
structure in wavelet helps the reconstruction and secondly, the 
Bayesian framework developed in | [4^ is very robust to noise; 
it infers noise from the measurements. 

B. RGB Images 

Next we consider the RGB images captured by the RGB 
sensor, and now the resolution is 217 x 302 x 3. The sensing 
matrix is constructed from rows of a Hadamard matrix of 
order N = 2^^ and the first 65534 elements are used. The 
scene is the real scene of four books as shown Q, | |20l . The 
reconstruction result is shown in Figure with diverse CSr. 

Note that by using compressive measurements, we can 
save the sensors as well as the bandwidth. As stated earlier, 
we may progressively get better results by receiving more 
measurements. One of the main usage of compressive sensing 
is to get features in limited data by using a small bandwidth. 
From the results in Figure we may identity high quality 
features from the reconstructed image at CSr around 0.1. If 
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Proposed 


CSr = 0.05 CSr = 0.1 CSr = 0.15 CSr = 0.2 CSr = 0.25 CSr = 0.3 










Proposed 


NLR-CS 


DAMP 


TVAL3 


GAP 


CSr = 0.05 CSr = 0.1 CSr = 0.15 CSr = 0.2 CSr = 0.25 CSr = 0.3 




Fig. 8. Real data: reconstruction results at different CSr with the diverse 
algorithms. The image is of size 128 x 128. Two photos (top: Lena, bottom: 
Alexander Graham Bell) are used as the scene. 


we want to get some details, for example, the book titles in 
Figure]^ we may need CSr around 0.2. On the other hand, if 
we only need to identify that these are “books” in Figure 
CSr at 0.05 may be sufficient. 

VII. Conclusions 

A novel compressive sensing reconstruction algorithm is 
developed via exploiting the low-rank property of overlap¬ 
ping patches. A general iteratively two-step framework for 
compressive sensing recovery is proposed. A denoising op¬ 
erator is used to update the estimate of the desired image 
(obtained by the projection of the measurements), which can 
be implemented by investigating the sparsity or low-rank 


CSr: 0.05 CSr: 0.07 CSr: 0.09 



Fig. 9. Real data: reconstruction results at different CSr with the proposed 
algorithm. The image is of size 217 x 302 x 3. 


property of the image patches. We develop a probabilistic 
regime by representing each patch via a Gaussian mixture 
model and impose low-rank on each Gaussian component to 
achieve the state-of-the-art compressive sensing reconstruction 
results, in particular when the measurement number is small. 
Additionally, the proposed low-rank GMM algorithm degrades 
to the low-rank piecewise linear estimator if each patch is 
modeled by a single Gaussian model. Extensive results on 
both simulation and real data demonstrate high performance 
of the proposed algorithm. 
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